% $Id: b_norm.m 320 2005-05-05 20:07:09Z odellj $
% b_norm - function version of boxy_norm.m - normalize continuous data from BOXY (unparsed) - function version
%
% elm 03/10/04 cleaned up print messages, supressed polyfit warnings
%               & chgd "filt' to 'norm'
% elm 2/24/11 added phase outlier cleaning
% elm 7/7/11 saving good_dist_chs, good_mrac_chs & good_mrdc_chs in output file

%cycle	exmux	A-AC	A-DC	A-Ph	B-AC	B-DC	B-Ph	time	group	mark	flag

function b_norm(good_dist_chs,log_file)
global hdr file_name i_block dc ac ph adc aac aph msg mod_freq id_tags
% temp actmp

b_norm_id = '$Id: b_norm.m 320 2005-05-05 20:07:09Z odellj $';
check = strcmp(b_norm_id,id_tags); % Make sure we don't already have the tag
if sum(check) == 0
    id_tags{size(id_tags,1)+1,1} = b_norm_id;  % Add it in then.
end

n_bad_points=12;    % number of anomalous initial data points to zap

in_path=[hdr.data_path hdr.exp_name 'opt\'];

[status,msg] = mkdir([hdr.data_path hdr.exp_name], ['norm00-00']);

out_path=[hdr.data_path hdr.exp_name 'norm00-00\'];

in_file_flg=fopen([in_path file_name]);
out_file_flg=fopen([out_path file_name]);
fclose('all');

% Open the log file if it exists
%
if ~strcmp(log_file,'')
    logid = fopen(log_file,'a');
else logid = 1; end

if ((in_file_flg>0) & (out_file_flg<0))
    
	start_row=84;
	stop_row=start_row+hdr.n_points*hdr.n_srcs-1;
	start_col=2;                            % Skip "cycle	exmux"
	stop_col=start_col+3*hdr.n_dets-1;        % don't bother reading "time	group	mark	flag"
	
	dc=ones(hdr.n_points,hdr.n_chans);
	ac=ones(hdr.n_points,hdr.n_chans);
	ph=ones(hdr.n_points,hdr.n_chans);
	tmp=ones(hdr.n_points,hdr.n_chans);
    
	%gain=load([in_path 'gains.txt']);          % get the gain settings for each detector

	% Get the Boxy data
	%
	fprintf(logid,'Reading %s\r\n',[in_path file_name]);
    
    line='';
    fid=fopen([in_path file_name]);
    line_num=0;
    while strncmp(line,'#DATA BEGINS',12) ~= 1;    
        line=fgets(fid);
        line_num=line_num+1;
    end
    
    line=fgets(fid); % read lables
    line=fgets(fid); % skip blank line
    line=fgets(fid); % read 1st data line
    % !!! 1st col is ch # !!!
    
    % figure out how many 'aux' channels there are
    %
    temp=sscanf(line,'%f');
    n_cols=length(temp); % length of 1st mux chan; includes aux data
    fclose(fid);

    start_row=line_num+2;
    stop_row=start_row+hdr.n_points*hdr.n_srcs-1;
    start_col=2;                        % Skip "cycle	exmux"
    % read extra columns "time,group,step,mark,flag,aux1-auxn,digaux"
    % which only exist on 1st mux channel
    % 'dlmread' will put zeros in missing cells !
    stop_col=start_col+3*hdr.n_dets-1;

    tmp_data = dlmread([in_path file_name],'\t',[start_row start_col stop_row n_cols]);
    data=tmp_data(:,1:3*hdr.n_dets); % optical data
    aux_data=tmp_data(1:hdr.n_srcs:end,stop_col+1:end);
    
    % copy desired extra columns to named variables
    %
    % time=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2)-7);
    % group=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2)-6);
    % step=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2)-5);
%     mark=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2)-4);
%     boxy_flags=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2)-3);
%     aux1=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2)-2);
%     aux2=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2)-1);
%     digaux=tmp_data(1:hdr.n_srcs:end,size(tmp_data,2));

	% Convert Boxy format to dc, ac & ph arrays
    %
    fprintf(logid,'Converting Boxy format the new way:\n');
    
    tmp=data(:,1:3:end); % extract ac data
    tmp=reshape(tmp,hdr.n_srcs,hdr.n_points,hdr.n_dets);
    tmp=permute(tmp,[2,1,3]);
    tmp=reshape(tmp,hdr.n_points,[]);
    ac(:,good_dist_chs)=tmp(:,good_dist_chs);
    
    tmp=data(:,2:3:end); % extract dc data
    tmp=reshape(tmp,hdr.n_srcs,hdr.n_points,hdr.n_dets);
    tmp=permute(tmp,[2,1,3]);
    tmp=reshape(tmp,hdr.n_points,[]);
    dc(:,good_dist_chs)=tmp(:,good_dist_chs);
    
    tmp=data(:,3:3:end); % extract ph data
    tmp=reshape(tmp,hdr.n_srcs,hdr.n_points,hdr.n_dets);
    tmp=permute(tmp,[2,1,3]);    
    tmp=reshape(tmp,hdr.n_points,[]);
    ph(:,good_dist_chs)=tmp(:,good_dist_chs);
    
    % Fix bad first points (Dennis' problem :)
    %
    for index = 1:n_bad_points
        dc(index,:)=dc(n_bad_points+1,:);
        ac(index,:)=ac(n_bad_points+1,:);
        ph(index,:)=ph(n_bad_points+1,:);
    end
    
	% fix phase wrap
    % 1st estimate mean phase of point 1 to 50
	% if a pont id more than 90 degrees different from the mean,
    %  add or subtract 360 degrees
    %
	fprintf(logid,'Fixing phase wrap\n');
	for i_chan=good_dist_chs
        if mean(ph(1:50,i_chan)) < 180
            wrapped_pts=find(ph(:,i_chan)>270);
            ph(wrapped_pts,i_chan)=ph(wrapped_pts,i_chan)-360;
        else
            wrapped_pts=find(ph(:,i_chan)>270);
            ph(wrapped_pts,i_chan)=ph(wrapped_pts,i_chan)+360;
        end
	end
    
    % Detrend (polynomial) phase data (fix Emily's problem ;-)
    % (also sets mean = 0 !)
    %
    fprintf(logid,'Detrending\n');
    y=1:hdr.n_points;
    x=y';
    warning off MATLAB:polyfit:RepeatedPointsOrRescale
    
    for i_ch=good_dist_chs
        poly_coeffs = polyfit(x,ph(:,i_ch),3);   % 3 => 3rd order
        tmp_ph = ph(:,i_ch)-polyval(poly_coeffs,x);
        ph(:,i_ch)=tmp_ph;
    end
    
	% Compute means
	%
	mrdc=mean(dc);
    dcneg=find(mrdc<0); % find chs w -dc
    mrdc(dcneg)=0; % and set them to 0
	mrac=mean(ac);
	mrph=mean(ph);
    
    % apply DC & AC quality thresholds
    %
    good_mrac_chs=find(mrac>=hdr.qual.norm_mrac_thr);
    bad_mrac_chs=find(mrac<hdr.qual.norm_mrac_thr);
    good_mrdc_chs=find(mrdc>=hdr.qual.norm_mrdc_thr);
    bad_mrdc_chs=find(mrdc<hdr.qual.norm_mrdc_thr);
    
    % set bad mrdc & mrac chans = 0
    %
    dc(:,bad_mrdc_chs)=0;
    ac(:,bad_mrdc_chs)=0;
    ph(:,bad_mrdc_chs)=0;
    dc(:,bad_mrac_chs)=0;
    ac(:,bad_mrac_chs)=0;
    ph(:,bad_mrac_chs)=0;
    
	% Normalize, and remove mean ('1') for filtering
	%
	fprintf(logid,'Rescaling (''Normalizing'')\n');
	for i_chan=good_dist_chs
        if mrdc(i_chan) ~= 0
			dc(:,i_chan)=dc(:,i_chan)/mrdc(i_chan)-1;
        end
        if mrac(i_chan) ~= 0
            ac(:,i_chan)=ac(:,i_chan)/mrac(i_chan)-1;
        end
        ph(:,i_chan)=(ph(:,i_chan)-mrph(i_chan));
	end
	ph=1e12*ph/(360*mod_freq); % convert phase to ps
    
    % remove phase (delay) outliers
    %
    SDph=std(ph(n_bad_points:end,:));
    ph_out_thr=hdr.qual.ph_out_thr; 
    for i_chan=good_dist_chs
        outliers=find(abs(ph(:,i_chan))>ph_out_thr*SDph(i_chan));
        if length(outliers)>0
            if outliers(1)==1;outliers=outliers(2:end);end % can't interp 1st pt
            if outliers(end)==hdr.n_points % can't interp last pt
                outliers=outliers(1:end-1);
            end
            n_ph_out(i_chan)=length(outliers);
            for i_pt=1:n_ph_out(i_chan)
                j_pt=outliers(i_pt);
                ph(j_pt,i_chan)=(ph(j_pt-1,i_chan)+ph(j_pt+1,i_chan))/2;
            end
        end
    end
    
    % save the cleaned up and normalized data
    % could save outliers, would need to make cell array of channels...
    %
    hdr_norm=hdr; % save hdr WITH UNIQUE NAME!
	msg=sprintf('Saving %s...\n',[out_path file_name]);
    fprintf(logid,'Saving %s\n',[out_path file_name]);
   	save([out_path file_name],'dc','ac','ph','mrdc','mrac','mrph','SDph',...
        'aux_data',...
        'poly_coeffs','id_tags','mod_freq','hdr_norm',...
        'good_dist_chs','good_mrac_chs','good_mrdc_chs','-mat');
% 'mark','boxy_flags','aux1','aux2','digaux'

else % if various files exist or don't...
    if (in_file_flg<0)
        msg=sprintf('%s does not exist.\n',[in_path file_name]);
        fprintf(logid,'%s does not exist.\n',[in_path file_name]);
    end
    if (out_file_flg>0)
        msg=sprintf('%s already exists!\n',[out_path file_name]);
        fprintf(logid,'%s already exists!\n',[out_path file_name]);
    end
end

if logid ~= 1; fclose(logid); end

% $Log: b_norm.m,v $
% Revision 1.26  2004/09/09 20:53:34  edmaclin
% bug 21 - replaced "'" with ";"  in"if logid ~= 1; fclose(logid)"
%
% Revision 1.24  2004/08/25 22:14:48  odellj
% Fixed the bad points code that wasn't running and also included the new id_tag structure code.
%
% Revision 1.23  2004/07/27 16:01:28  edmaclin
% zap "n_bad_points" initial points
%
% Revision 1.22  2004/07/15 17:39:16  edmaclin
% fclose 'all', not specific fids (which may not exist)
%
% Revision 1.21  2004/07/13 22:26:47  odellj
% Added the new variable _id to the subprogram to output to the inidiviually created files to identify the version of the producing function.
%
% Revision 1.20  2004/07/12 20:00:53  edmaclin
% *** empty log message ***
%
% Revision 1.19  2004/07/12 19:59:08  edmaclin
% saves poly_coeffs in output file
%
% Revision 1.18  2004/07/09 21:51:31  odellj
% Updated the logging of the subroutines.  Added log_file as a parameter for each of their calls.
%
% Revision 1.17  2004/07/08 19:14:31  edmaclin
% fixed typo in bad chans code - elm
%
% Revision 1.16  2004/07/07 18:58:10  edmaclin
% fixed non-zeroed bad chans - elm
%
% Revision 1.15  2004/07/07 18:07:42  odellj
% Prepped for logging.
%
% Revision 1.14  2004/06/30 18:13:41  odellj
% Added more verbose fprintf statements.
%
% Revision 1.13  2004/05/19 21:46:25  edmaclin
% good_chs not good_chans
%
% Revision 1.12  2004/05/19 20:36:10  edmaclin
% use .gch not .dis
%
% Revision 1.11  2004/05/13 17:49:16  edmaclin
% returns pico not nano seconds for phase
%
% Revision 1.10  2004/05/11 19:50:20  edmaclin
% convert phase to pico seconds (mod_freq is global)
%
% Revision 1.9  2004/05/11 17:52:10  edmaclin
% convert phase to pico seconds (mod_freq is global)
%
% Revision 1.8  2004/05/07 22:16:16  edmaclin
% declared SD arrays to avoid missing (n=0) bins
%
% Revision 1.7  2004/03/31 22:06:40  edmaclin
% changed filt to norm
%
% Revision 1.6  2004/03/29 22:11:56  odellj
% Changed norm to filt.
%
% Revision 1.5  2004/03/29 16:53:59  odellj
% Fixed the '\' issue.
%
% Revision 1.4  2004/03/23 22:35:40  odellj
% Removed old way of converting Boxy format.
%
% Revision 1.3  2004/03/16 20:47:25  edmaclin
% vectorized conversion
%
% Revision 1.2  2004/02/05 20:40:51  odellj
% Added Id and Log tags.
